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Abstract 

y(t,p] 

This paper describes and illustrates two matched -filter- 
theory based schemes for obtaining maximized and time 
correlated gust-loads for a nonlinear airplane. Hie first scheme y^(t) 
is computationally fast because it uses a simple one-dimensional 
search procedure to obtain its answers. 1 he second scheme is yy(t) 
computationally slow because it uses a more complex multi- 
dimensional search procedure to obtain its answers* but it 
consistently provides slightly higher maximum loads than the 7 (t) 
first scheme. Both schemes are illustrated with numerical 
examples involving a nonlinear control system. z y (t) 


time response of output quantity y to iterative 
excitation waveform 

time response of output quantity y to an impulse 

time response of output quantity y to excitation 
matched to y 

time response of output quantity z 

time response of output quantity z to excitation 
matched to y 


Nomenclature 


scalar multiplying factor used in the KS function 


BM 

bending moment 

CXC(l,p) 

iterative excitation waveform time history 
from optimization design variables (p) 

F(p) 

objective function 

g(p) 

constraint function 

k 

impulse strength 

KS(p) 

Kreisselmeier Steinhauscr function 

i 

iteration number 

P 

vector of design variables 

RMS 

root mean square 

RMS(x) 

root mean square of quantity x 

t 

iinx k 


torsion moment 

lo 

arbitrary time shill 

WOBTM 

wing outboard torsion moment 

WRBM 

wing root bending moment 

yd) 

time response of output quantity y 
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Og root-mean -square of gust velocity 


Introduction 

Time-correlated gust loads are time histories of two or 
more load quantities due to the same disturbance time history. 
Time correlation provides knowledge of the value (magnitude 
and sign) of one load when another is maximum. Figure 1 
contains an illustration of time -correlated gust loads. At least 
two analysis methods have been identified (refs. 1 and 2) that are 
capable of computing maximized time-correlated gust loads for 
linear aircraft. Both methods solve for the gust profile (gust 
velocity as a function of time) that produces the maximum load at 
a given location on a linear airplane. Time-correlated gust loads 
are obtained by re-applying this gust profile to the airplane and 
computing multiple simultaneous load responses. The load 
responses are physically realizable at all instants of time during 
the time history (including the instant at which a given load is 
maximum), and correlated loads at any of these instants may be 
applied to aircraft structures. 



Figure 1.- Time -correlated gust loads. 
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Within the past several years there has been much interest 
in obtaining a practical analysis method that is capable of solving 
the analogous problem for nonlinear aircraft. Such an analysis 
method, and the broader issues of harmonizing existing criteria 
and minimizing the number of different analyses required for 
certification, have been the on going focus of an international 
committee of gust loads specialists formed by the U S Federal 
Aviation Administration (ref 3). 'I Tie method ol reference I is 
capable of being, but so far has not been, applied to nonlinear 
aircraft. As stated in reference I, to make the method practical 
for nonlinear applications, a search prtxrcdurc is essential. The 
mctluxl of reference 2 is based on Matched Filter Ilieory (MIT) 
and, in its current form, is applicable to linear systems only. 

The purpose of the present p;i|x.T is to extend the mctluxl 
described in reference 2 to nonlinear systems. The extension 
uses MFT as a starting point and then employs a constrained 
optimization algorithm to attack the nonlinear problem. 
Reference 4 contains a status report (as of April 1990) on the 
development of this extension. In the present paper, 
methodology development is refined and completed and 
numerical results are presented that illustrate an application of the 
method to an airplane with a nonlinear control system 

Review of Matched Filter Theory For Linear Systems 

As background for the remainder of this paper, this 
section presents the key results from reference 2, the original 
development of MFT for computing maximized and lime- 
correlated gust loads for linear airplanes. 


Figure 2 depicts the steps that must be employed to 
implement MFT and illustrates the intermediate and final 
products of the process. Frequency-response-function 
representations of atmospheric turbulence and airplane loads are 
combined in series and represent the "known dynamics" boxes 
in the figure. One-dimensional Gaussian atmospheric 
turbulence with the von Karman power spectrum is chosen. 
I x>ad y is the load to be maximized. Loads z\ through Zn are the 
loads to be time correlated with load y. There are three major 
steps in the process: 

Step 1 The application of an impulse function of 
unit strength to the combined linear system, 
producing the impulse response of load y. 

Step 2 The normalization of this impulse response 
by its own RMS value, followed by its 
reversal in time. 

Step 3 The application of this normalized reversed 
signal to the combined linear system, 
producing time histories of load y and time 
histories of loads Z| through z n . Within the 
time history of load y, the maximum value is 
ymax- 

For simplicity of discussion throughout this paper, these three 
steps will be referred to as "linear MFT." Within each step, 
significant results, interpretations, and implications are: 
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Within Step 1 The impulse response of atmospheric 
turbulence is shown as an intermediate 
quantity. The impulse responses of 
loads z\ through can be computed 
but are not used. 

Within Step 2 The normalized reversed signal, named 
the matched excitation waveform, now 
has a unit RMS value. MIT guarantees 
that there is no other unit RMS 
excitation waveform that, when applied 
to the combined linear system, will 
produce a v\..uc ol y larger than y n \ax 
This guarantee is a fundamental result 
of Ml T 

Within S tep 3 The response of atmospheric turbulence 
to the matched excitation waveform, 
named the critical gust profile, has a 
given level of probability. The 
responses of atmospheric turbulence to 
all other unit RMS excitation 
waveforms have the same given level of 
probability. Load y is named the 
matched load and loads zj through z n 
arc time correlated with load y. At each 
instant ol time, the time history of the 
matched load is proportional to the 
autocorrelation function of load y from 
a conventional random process 
analysis. At each instant of time, the 
time histories of loads /.| through are 
proportional to their corresponding 
cross-correlation functions with load y. 

Tor more detailed information about the application of MIT to 
the calculation of maximized and time -correlated gust loads, 
including background, theoretical development, numerical 
implementation, and example calculations, the reader is urged to 
consult reference 2. 


The Application of Matched Filter Theory to 
Nonlinear Systems 

The goal of MhT as applied to nonlinear systems is the 
same as that for linear systems: to find the maximized response 
time history, the maximum value of the response within that lime 
history, and the time correlated response time histories Because 
the system is now nonlinear, superposition and other convenient 
features of linear MIT no longer hold, the most convenient 
being that problems can be solved directly. l ; or nonlinear 
systems, problems cannot be solved directly and a search 
procedure is the only practical means of finding the solution. 
Here ’’the solution” means searching for and finding the 
excitation waveform that maximizes y max . The search is 
conducted systematically, subject to the constraint that the 
excitation waveform have unit RMS, using the techniques of 
numerical optimization. 

Because superposition no longer holds, the magnitude 
and character of the responses are not necessarily proportional to 
the magnitude of the input. I or the remainder of this paper two 
input magnitudes are important: k, the strength of the initial 
impulse; and Gj,, the design value of the gust intensity. (Lor 
linear Mf ; T, the magnitude of both of these quantities was 
unity.) l ; or nonlinear systems the shape of the excitation 
waveform is a function of k and a p , and, consequently, the 
quantity y n iax is also a function of these parameters. 

This section of the paper is divided into two subsections. 
The first describes a one-dimensional search procedure and the 
second describes a multi dimensional search procedure. The 
original goal of the present research was to develop the multi- 


dimensional search only. However, in an effort to choose a 
reasonable starting point for the multi -dimensional search the 
authors happened upon the one-dimensional search, and as will 
be shown in the numerical results section of this paper, the one- 
dimensional search turned out to be a valuable procedure in its 
own right. 

One-Dimensional Search Procedure 

The one-dimensional search procedure performs a 
systematic variation of the quantity k to find the shape of an 
excitation waveform that maximizes y ma x. Figure 3 contains a 
signal flow diagram for this search procedure. Figure 3 is very 
similar to figure 2, but contains some subtle yet important 
differences that arc indicated by the shaded boxes and by 
quotation marks. In figure 3 the initial impulse has a non-unity 
strength; the aircraft loads portion of the known dynamics box 
contains nonlinearities; and the shape of the excitation waveform 
and the value of y ma x are functions of the initial impulse 
strength. In addition, the "matched” excitation waveform and 
the "matched" load are shown in quotes because, for nonlinear 
systems, there is no guarantee that y^ is a global maximum. 

The application of the one-dimensional search procedure 
is as follows: 


Step 1 Select a design value of Gg. 

Step 2 Select a range of values of k. 

Step 3 For each value of k, perform steps 1 
through 3 of linear MFT, obtaining values 
of ymax and corresponding "matched" 
excitation waveforms. 

Step 4 From step 3 above, find the maximum 
value of y m ax and its corresponding 
"matched” excitation waveform. 

Mulli-Dimensional Search Procedure 

The multi dimensional search procedure uses as its 
starting point the maximum value of y max and the corresponding 
"matched" excitation waveform obtained by performing the one- 
diinensional search procedure. In an attempt to obtain an even 
larger value of y max > the shape of the "matched" excitation 
waveform is further optimized subject to the unit-RMS constraint 
mentioned above. The major features of this optimization 
procedure arc illustrated in figure 4. 

Initial excitation waveform. - Beginning in the upper left 
corner ol figure 4, the box labelled initial excitation waveform is 
the waveform resulting from the one-dimensional search 
procedure. 

Coefficient generator. - The purpose of the "coefficient 
generator" is to approximate the initial excitation waveform by a 
finite number of coefficients. These coefficients are first used as 
the design variables in the optimization procedure and later used 
to "reassemble" the final optimized waveform time history. 
Because the execution lime of the optimization procedure is 
roughly proportional to the number of design variables, an 
important consideration was to choose an approximation 
technique that produced in the smallest number of coefficients 
for a given goodness of approximation. Several approximation 
techniques were considered including the use of orthogonal 
functions and cubic splines. 

The orthogonal functions considered were Fourier series 
and Chebyshev polynomials. It was found (ref. 4) that it took 
about an order of magnitude fewer coefficients to approximate 
the initial excitation waveform with Chebyshev polynomials than 
to approximate it with Fourier series. Figure 5 shows an initial 
excitation waveform and approximations to it using first a 20- 
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Figure 3.- Schematic for one-dimensional search procedure. 



Figure 4.- Schematic for multi -dimensional search procedure. 
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term and then a 54- term Chebyshev approximation. The choice 
of how many coefficients to include in the approximation was 
made by examining a plot of the absolute value of the coefficient 
versus the coefficient number. 


Normalized 
Excitation I— 
Waveform I 

0 


Normalized 
Excitation 
Waveform 

0 5 Time, sec. 10 15 

(b) Chebyshev approximation using 20 terms. 


Normalized 
Excitation 
Waveform 

0 5 Time, sec. 10 15 

(c) Chebyshev approximation using 54 terms. 

Figure 5.- Excitation waveform approximations using 
Chebyshev polynomials. 




decreases when the objective function increases and increases 
when the constraint is violated. Here, the optimizer varies the 
coefficients to minimize the KS function. 


KS(p) = ^ log e [exp(p { (gfa. i )+ 1 - } ) 

+ ex P(P g(Pi))] (3) 

where F and g are defined in equations (1) and (2). 

The shaded box in figure 4 contains the detailed 
sequential analysis steps performed. The optimizer cycles 
through the shaded box many times per iteration to obtain 
gradient information. The accompanying illustrations in the 
figure (boxes with rounded comers) show intermediate results 
and the nonlinear model. 

Nonlinear simulation. The most computationally 
expensive aspect of this multi-dimensional search is the 
nonlinear simulation, shown in the shaded box in figure 4. The 
original simulation model was written in MATRIX* 
SYSTEM BUILD simulation language (ref. 7). This high-level 
language proved very slow considering the large number of 
simulations required by the optimization algorithms. Therefore, 
the entire process was converted to FORTRAN with the code for 
the simulation model being created by MATRIX* 
HYPER BUILD (ref. 8). This greatly increased throughput so 
that optimization results could be obtained in a reasonable 
amount of time. In addition, double precision was used so that 
sufficiently accurate gradient information could be obtained from 
the simulations. 


A waveform approximation technique based on cubic 
splines was also considered. For this method, discrete points in 
the waveform time-history were selected as design variables. 
The amplitudes of these points could then be changed by the 
optimization routine. By using a cubic spline the waveform time 
history can be assembled from these discrete points as needed in 
the optimization algorithm. As will be shown, this technique did 
not provide as accurate an approximation as Chebyshev 
polynomials, yet the number of design variables could be 
substantially reduced. 

C onstrained optimization algorithms. - Two optimization 
algorithms were considered: a direct method and an indirect 
method. This section briefly describes these algorithms. 

The direct method (ref. 5) chooses search directions 
(systematic variations of the values of the coefficients) that 
minimize an objective function without violating constraints. 
The objective of this particular application is to maximize one of 
the loads. This objective function is formulated such that 
increasing the load, maximizing wing root bending moment 
(WRBM) for instance, reduces the objective function. 

Objective: minimize F(p) = - { max[ y(t,p) | }2 (!) 

over all t 

The only constraint is that of keeping the RMS of the 
excitation waveform, exc(t,p), at or below a specified value. In 
actuality the waveform RMS should remain constant, implying 
the necessity for an equality constraint. However, in practice an 
inequality constraint performs very well, with the optimizer 
driving the value of the constraint to its maximum-allowable 
value. This, in effect, gives the desired result while not overly 
constraining the design space. 

Constraint: g(p) = {RMS| exc(t,p) | - 1.0} ^0.0 (2) 

The indirect method (ref. 6) converts a constrained 
optimization problem into an unconstrained problem by 
combining the constraints with the objective functions into a 
single composite function known as the Kreisselmeier- 
Steinhauser function. The KS function is formulated such that it 


Numerical Results 

This section presents numerical results for the one- 
dimensional and multi-dimensional search procedures. 

Example System 

A nonlinear simulation model of the ARW-2 drone 
aircraft equipped with a nonlinear control system was 
constructed to demonstrate both search procedures. Figure 6 
shows a block diagram of the simulation model and includes the 
aeroelastic plant, a gust load alleviation (GLA) control law, and 
nonlinear control elements. The aeroelastic plant is a linear, s- 
plane aeroelastic half-model consisting of two symmetric rigid- 
body modes and three flexible modes. Unsteady aerodynamics 
were obtained using the doublet lattice method (ref. 9). The 
model also includes the dynamics of the control surface 
actuators. The two-input/two-output GLA control law was 



Figure 6.- Nonlinear model of aircraft and gust load alleviation 
system. 
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obtained using a Linear Quadratic Gaussian design approach 
with the intent of reducing wing root bending moment (ref. 10). 
The nonlinear elements impose deflection limits of ±KT on the 
elevators and 0 to 4 10° on the ailerons, simulating spoilers. 

Two versions of this model will be considered: a simple 
version having only rigid-body modes (10 states) that served as 
a computationally economical test case; and a much larger, fully 
flexible model (32 states). These models will be hereinafter 
referred to as the rigid and flexible models, respectively. 
Analysis conditions for this model are 0.86 Mach number and 
24,000 feet altitude. 


One-Dimensional Search 

The one -dimensional search procedure was performed on 
the rigid model. In order for the control-surface deflection limits 
to become active in the simulation, an artificially large value of 
gust intensity, cjg, was required. With this value chosen, values 
of the initial impulse strength, k, were selected over the range 
from 0.1 to 1000. Typical results, in this case corresponding to 
maximizing wing root bending moment (WRBM), are presented 
in figure 7. Figure 7 contains 3 plots. The lower portion of the 
figure contains a plot of WRBM as a function of initial impulse 
strength: the solid circles are actual numerical results; the solid 
line is a faired line. Three sets of time histories of excitation 
waveform and its associated maximized WRBM are shown in 
the upper portion of figure 7. The corresponding values of k are 
indicated on the figure. For this particular example, a value of k 
of about three yields the largest maximized value of WRBM. 



Figure 7 - Results from one -dimensional search procedure. Rigid model. 


Multi-Dimensional Search 

Maximizing WRBM in the rigid model will again be 
considered. To more clearly demonstrate the method, a 
nonoptimal value of k (ki) was chosen to obtain the initial 
excitation waveform for this example. The Chebyshev 
waveform approximation was utilized in this example and 
required 54 coefficients (design variables). The direct optimizer 
was used and the resulting convergence history is shown in 
figure 8. This figure shows that the optimization scheme 
succeeds in increasing the maximum WRBM while ultimately 
satisfying the RMS constraint on the excitation waveform. 

Figure 9 shows a comparison of the initial and optimized 
time histories for this test case. By considering figures 8 and 9, 
two conclusions can be drawn: (1 ) the multi-dimensional search 
procedure succeeds in altering the shape of the excitation 
waveform which causes alterations of the various time-histories 
and an increase in the maximum WRBM; and (2) the procedure 
achieves this increase in maximum load while maintaining the 
RMS constraint on the excitation waveform. 


Figure 10 presents the results of applying the multi- 
dimensional search procedure a number of times. Again WRBM 
is the load of interest and again the rigid model is employed. 
The direct and indirect optimizers were exercised and excitation 
waveforms were approximated by Chebyshev polynomials and 
by cubic splines. The numbers in parentheses indicate the 
number of terms required to obtain an acceptable approximation 
to the waveform. The solid line in figure 10 is the faired line 
presented originally in figure 7. For the example shown, the 
direct optimizer with the Chebyshev polynomials was the best 
combination for reproducing the results from the one- 
dimensional search procedure and for obtaining the largest 
maximized loads. 

The squares and circles in figure 10 represent results of 
the initial and final values of the maximum WRBM. In all cases, 
at each value of k, the lower symbols represent the initial value 
of maximized WRBM and the upper symbols represent the final 
value. Regardless of the starting point, the direct optimizer 
consistently achieved the same value of maximized WRBM, 
indicating that it was more insensitive to the starting point than 
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was the indirect optimizer. The direct optimizer also yielded 
larger maximum load values and converged in fewer iterations. 
In addition, the direct optimizer was able to converge with both 
Che by she v and cubic spline approximations to the excitation 
waveform while the indirect optimizer was not able to converge 
when the spline approximation was used. 

In figure 10 the lower gray squares are below the solid 
line but the upper gray squares are very close to the highest value 
of WRBM achieved. The lower symbols being below the solid 
line is due to the fact that the spline approximation does a poorer 
job at reproducing the initial excitation waveform than does the 
Che by she v approximation. However, the spline approximation 
does a very adequate job of reproducing the final excitation 
waveform and produces maximized WRBMs very close to those 
produced by the Chebyshev approximation. 

From the data presented in figure 10, the usefulness and 
value of the one-dimensional search procedure becomes 
obvious. A very simple and straightforward variation in a single 
quantity, k, produced a value of maximized load which is within 
1.1% of the value produced by the significantly more complex 
multi-dimensional search procedure. 
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Figure 9- Comparison of initial and optimized time-histories. Rigid model. 
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Figure 10 - Summary of results for maximizing WRBM using 
the one-dimensional and several multi dimensional 
searches. Rigid model. 



Figure 1 2.- Summary of results for maximizing WOBTM 
using the one-dimensional and several multi- 
dimensional searches. Flexible model. 


Figure 11 shows results for maximizing WRBM using 
the flexible model. The solid line again represents results from 
the one-dimensional search procedure and the squares represent 
results from the multi dimensional search. In this example, 
because of the higher-frequency dynamics in the excitation 
waveform, 124 Chebyshev polynomials were required to 
produce an adequate approximation of the waveform. The 
indirect optimizer was not considered here for the reasons cited 
above. Also, the spline approximation was unsuccessful for this 
model, even with the direct optimizer. Consistent results were 
obtained for the flexible model only when the direct optimizer 
and the Chebyshev polynomial approximation were used. 

The usefulness and value of the one-dimensional search 
procedure is again illustrated by Ihe results from figure 1 1. The 
simple procedure produced a value of maximized load which is 
within 0.78% of the value produced by the multi dimensional 
search procedure. 

So far, the only load considered has been wing root 
bending moment. Figure 1 2 shows results for maximizing wing 
outboard torsion moment (WOBTM) for the flexible model. The 
solid line represents results from applying the one-dimensional 
search procedure to WOBTM. The general shapes of the solid 
lines in figures 10, II, and 12 are very similar. For the multi- 
dimensional search procedure, 220 Chebyshev polynomials 



Initial Impulse Strength, k 

Figure 1 1 Summary of results for maximizing WRBM using 
the one -dimensional and several multidimensional 
searches. Flexible model. 


were required to produce an adequate approximation of the 
excitation waveform. The solid squares represent only the initial 
values of maximized WOBTM; final values are not shown 
because, due to the large number of coefficients, a converged 
solution was not obtained. 

Based on the ability (demonstrated in figures 10 and 11) 
of the one-dimensional search procedure to produce maximized 
loads very close to those produced by the multi-dimensional 
search procedure, it is the opinion of the authors that the value of 
the maximum WOBTM is very close to the value predicted by 
the one-dimensional search. 


Concluding Remarks 

This paper has presented two search procedures (one- 
dimensional and multi-dimensional) for obtaining maximized 
gust loads for a nonlinear airplane using Matched Filter Theory 
and constrained optimization. The procedures were applied to an 
example airplane with a nonlinear gust load alleviation active 
control system. Numerical results confirm that the procedures 
are successful in maximizing the loads of interest. For the 
example configuration considered, the one-dimensional search 
succeeded in obtaining loads that were within about 1% of the 
loads produced by the multi-dimensional search procedure. This 
suggests that the simpler one-dimensional search method may be 
the most practical implementation. In addition, further study is 
required of the modelling of excitation waveforms characterized 
by high-frequency contents. 
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